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Abstract 

We present in this paper a multigroup model for radiation hydrodynamics 
to account for variations of the gas opacity as a function of frequency. The 
entropy closure model (Mi) is applied to multigroup radiation transfer in a 
radiation hydrodynamics code. In difference from the previous grey model, 
we are able to reproduce the crucial effects of frequency- variable gas opacities, 
a situation omnipresent in physics and astrophysics. We also account for the 
energy exchange between neighbouring groups which is important in flows 
with strong velocity divergence. These terms were computed using a finite 
volume method in the frequency domain. The radiative transfer aspect of the 
method was first tested separately for global consistency (reversion to grey 
model) and against a well established kinetic model through Marshak wave 
tests with frequency dependent opacities. Very good agreement between the 
multigroup Mi and kinetic models was observed in all tests. The successful 
coupling of the multigroup radiative transfer to the hydrodynamics was then 
confirmed through a second series of tests. Finally, the model was linked to 
a database of opacities for a Xe gas in order to simulate realistic multigroup 
radiative shocks in Xe. The differences with the previous grey models are 
discussed. 
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1. Introduction 

The study of radiative transfer and its interaction with matter has an 
extremely wide range of applications ranging from medical imagery to astro- 
physics. In many cases, as for example in stellar atmospheres, the radiation is 
considered as a physical probe which provides access to the thermodynamical 
properties of the flow through the spectrum of emission and absorption lines. 
However, the radiation often has a very important dynamical role in the sys- 
tem. It cannot only be considered as a passive probe, but as an integral part 
of the equations governing the system dynamics. 

The equation of radiative transfer (ignoring scattering) is 



( c ^ + • '^^ H^, n, u) = (^5(x, t, u) - /(x, t; n, u)^ 



(1) 



where / is the specific intensity of the radiation, v the frequency, c is the 
speed of light, the absorption/emission coefficient and B the black body 
specific intensity, n, x, and t are the angular, spatial and temporal variables, 
respectively. As the radiation intensity depends on seven variables in three- 
dimensions, solving the full transfer equation coupled to the hydrodynamics 
to tackle radiation hydrodynamics (RHD) problems is still out of reach of 
modern computational architectures, even with the remarkable and constant 
increase in computing power. 

In order to overcome this difficulty, much effort has been spent in recent 
years developing mathematically less complicated, yet accurate approxima- 
tions to the equations of radiative transfer. Such approximations include 
diffusion approximations and moment models jsi [ill, [13, [2l| . All of 



these approximations use frequency and/or angle-integrated variables which 
greatly simplify the calculations. The approximations due to the angular 



integration have been widely studied (see for example Olson et al. 181]). 
However, in many situations, the quantities involved in the equations of 
radiative transfer (in particular the absorption and scattering coefficients) 
depend strongly on frequency, and the so called 'grey' approximation (in- 
tegrated over all frequencies) is no longer appropriate. Only very recently 
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have models which take into account variations in frequency been developed 
22I, M, Q, |20 |. The common practise is to split the frequency domain into 



a finite number of bins or groups and the equations of radiative transfer are 
solved within each group; this is known as a multigroup method. Such a 
scheme is then capable of allowing for gas opacity variations in the frequency 
domain providing a more accurate description of radiative transfer processes. 
The model we present in this paper is an extension of the moment model in- 



troduced in Turpault [22| which couples the frequency-dependent radiation 
to the hydrodynamics. 

Moment models are obtained by computing successive angular moments 
of the radiative transfer equation. One obtains a hierarchy of equations 
for moments of the specific intensity. Basically, each equation describes the 
evolution of the n*'^ moment as a function of the divergence of the (n-(-l)*'^ 
moment. For instance, the equations giving the evolution of the first two 
moments are 

dtE, + V-F, = a,{4nB-cE,) 

where Ey^Y^, and are, respectively, the radiative energy density, the ra- 
diative energy fiux, and the radiative pressure, which are defined in terms of 
the zeroth, first and second moments of the specific intensity as 

Eu = - (£ /(x, t; n, z/) dfl 

" 7 

Yy= i) n /(x, t; n, u) dVl (3) 

Pj, = - y n n /(x, t; n, z/) (in . 

The transfer equation is formally equivalent to an infinite hierarchy of 
moment equations. In order to have a tractable moment model, one must 
cut this hierarchy at some given order. A closure relation is then needed in 
order to express the moment of highest order as a function of the others. 

In this paper, we develop for the first time the coupling of the Mi moment 
model for radiative transfer to the hydrodynamics to create a multigroup 
model for RHD. A finite volume method used to compute the additional 
terms due to frequency variations is described. We then present a series of 
tests for the multigroup model studying first the radiative transfer alone and 
then the radiation coupled to the hydrodynamics. The strengths and future 
developments of the model are finally discussed. 
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2. The multigroup model for radiation hydrodynamics 



2.1. The monochromatic equations of radiation hydrodynamics 

The equations of radiation hydrodynamics describe the effects of radiative 
transfer on a moving fluid. The fluid evolution is determined by the classical 
conservation equations (mass, momentum, and energy) which are coupled to 
the radiative transfer equations (|2]) through source terms characterizing the 
momentum and energy exchanges between the fluid and the radiation. 

In order to write the RHD equations, one has to choose the frame in 
which to evaluate the radiative quantities: laboratory frame or comoving 
frame (i.e. the frame moving with the fluid). The laboratory frame is con- 
venient because the the left-hand side of the system remains hyperbolic and 



thus globally conservative [IJ] . However in this frame, interactions with mat- 
ter become complex because of Doppler and aberration effects that have to 
be incorporated in the source terms. On the other hand, using the radiative 
quantities expressed in the comoving frame 13|] adds non-conservative terms 
to the equations, and conversions of the radiative quantities between comov- 
ing and lab frames are required in order to be compared to observations as 
any measurement will almost certainly be carried out in the lab frame. How- 
ever the source terms coupling matter and radiation remain unaffected by 
the fluid motions. 

We have chosen to express radiative quantities in the comoving frame for 
the greater simplicity of the source terms. The equations of non-relativistic 



RHD (to order u/c) can then be written as |l3l . Il5l . |2| 



dtp 


+ V 


(pu) 


dt{pu) 


+ V 


(pu (g) u + pi) = 


dte 


+ V 


(u(e + p)) 







{a^/c)F^diy 

[ (a,{AnB-cE, 
(o-,,/c)u ■ du 



(4) 



dtE^ V ■ F^ + : Vu + V ■ (nE^) - 9^(z/P^) : Vu = a^(47r5 - cE^) 
dtF^ + c^V ■ + F^ ■ Vu + V ■ (u ® F^) - d^iuQ^) : Vu = -a^cF^ 

(5) 

where p is the gas density, u the velocity, e the total gas energy, p the gas 
pressure, and Qu is the third moment of the speciflc intensity 



n (g) n (g) n /(x, t; n, u) dQ 



(6) 
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The tensorial contractions are defined by P : Vu = Fijd^u^ and Q : Vu = 



2.2. The multigroup equations of radiation hydrodynamics 

Equations (jl]) and ([5]) are all Eulerian, but the radiative quantities are 
evaluated in the frame comoving with the fluid. In a grey model, system (jS]) 
is integrated from to oo in frequency and the terms involving frequency 
derivatives dv{v¥y) and dp^uQ^v) vanish. However, in a multigroup model 
these terms remain and are in fact of great importance; they govern energy 
transfers between neighbouring groups. 

In a multigroup model, the frequency domain is divided into a finite 
number of bins or groups and the radiative transfer equations are integrated 
and solved within each group. The integrals in the source terms of the 
hydrodynamic equations (jlj) then become sums of source terms over the total 
number of groups. Systems (jl]) and ([5]) become 



dtp 



+ V-(pu) 





Ng 



dt{pu) + V-(pu®u + pI) 




Ng 



(7) 



dte 




9=1 

{aFg/c)u ■ Fg 




= c 



{apgQg{T) - CTEgEg) 




(JFgCFg (8) 



with 




(9) 
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where X = E, F, ¥, Q which represent the radiative energy, flux, pressure 
and heat flux inside each group g which holds frequencies between i'g-1/2 
and Vg+i/2- Ng is the total number of groups and 6g(T) is the energy of the 
photons having a Planck distribution at temperature T inside a given group. 
The absorption coefficients apg, aEg and apg are the means of inside a 
given group weighted by the Planck function, the radiative energy and the 
radiative flux respectively. 

In order to integrate the previous system, it is necessary to introduce a 
closure relation giving P and Q as a function of E and F. The closure we 
have chosen is based on the Mi model and is presented below. 

2.3. The multigroup Mi model 

The Ml model |3] uses the first two moment equations to approximate 
the equation of radiative transfer. It has the great advantage over flux- 



limited diffusion models [16|, [121, IJJJ of being valid in both the diffusion and 
free-streaming limits while maintaining a directionality in the propagation 
of the radiation. Shadows can be created with the Mi method while flux- 
limited diffusion considers the radiative flux to always be colinear to the 
radiative temperature gradient, which can result in radiation propagating 
around corners js], 0] . 

In the Ml model the radiative pressure is expressed as P = 3E where 
D is known as the Eddington tensor. The expression for D is obtained by 
minimizing the radiative entropy which yields 

1 - Y 3y - 1 F ® F 
B= —^1 + ^ (10) 



where 

X = , (11) 

5 + 2v/4 - 3/2 

and / = is the ratio of the grey flux to the flux free-streaming limit, 
also known as the reduced flux. The quantities without the subscript u 
represent quantities integrated over the entire frequency range. Note that by 
definition of E and F, we have / < 1, which implies that the radiative energy 
is transported at most at the speed of light. In one dimension we simply have 
P = xE. We have plotted x iu Fig. [T]as a function of / (red). This closure 
relation recovers the two asymptotic regimes of radiative transfer. In the 
free-streaming limit (i.e. transparent media), we have / = 1 and x = 1- Ou 
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the other hand, in the diffusion hmit, / = and x = 1/3, which corresponds 
to an isotropic radiation pressure. 

In order to express the radiative heat flux Q as a function of the lower 
angular moments using the Mi closure, we define Q = MEc. Due to the sym- 
metry of the specific intensity distribution function around the axis defined 
by the direction of propagation of the radiative flux, it can be shown that 

EI = ipi{U6jk + ijSik + ikSij) + ip2(fifjfk) (12) 
in which fj are the components of the reduced flux vector f 



^1 



(/-2 + a)(/ + 2-a) 
4/(a-2)5 



12 In 



/-2 + a 

f + 2-a 



(f^ + 2af-7f-Aa + : 



+ A8f - 9a f - 80/ + 40a/ 



(13) 



and 



/3(a - 2)5 



60 In 



/-2 



/ 



16a + 32) + 5Aaf - 465/^ - 61^af + 2140/^ + 1056a/ - 2112/ 



(14) 



where a = a/4 — 3/^. We plot ipi and as a function of / in Fig. [T] (green 
and blue). Note that in one dimension, f = (/, 0, 0) and we have Q = ipEc 
where 

ilj = 3ipj + if2f (15) 

and ip is plotted in Fig. [1] (orange). This fully defines the evolution of the 
radiative energy and pressure of the model coupled to the hydrodynamics of 
the system. 

A natural way to extend this closure to a multigroup model would be to 
minimize the total radiative entropy, which is a rather complex procedure. 
However, Turpault j22| has shown that applying inside each group a closure 
formally equivalent to the Mi closure leads to almost indistinguishable re- 
sults; a strategy which we have therefore adopted for its greater simplicity. 
We define for each group the radiative pressure as Pg 



^gEg where 



IF l|2 



(16) 
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Figure 1: x (red), ipi (green), ip2 (blue) and -0 (orange) as a function of /. x, ipi and 
ip2 are symmetric with respect to the ordinates axis, ip is symmetric with respect to the 
origin. 



^9 = . , o /^A^ (17) 



(■2 

5 + 2^/4^ 

IIF II 

and fg = The heat flux Qg = MgEgC is computed in the same manner. 

2.4- A finite volume method for the frequency derivatives 

The only terms which were not included in our previous grey RHD models 
[3] are the terms in (jSj) involving the frequency differentials. In order to 
evaluate these terms, we adopt a finite volume method in the frequency 
dimension. We present here this method in the one-dimensional case, but 
its extension to several dimensions is trivial. Retaining only the time and 
frequency derivatives of the radiative energy and flux equations of system 
we obtain 

dtE,-V d,{uF,) = 

dtF,-Vd,{uQ,) = ^'""^ 

where P = V ■ u. We now assume that the frequency group boundaries 
(z/g-i-1/2) are equivalent to the volume elements' boundaries in the frequency 
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dimension. The finite volume discretization of f llSp gives 





(19) 



' ' - P ('^.+1/2Q;Vi/2 - ^.-l/2Q^l/2 j = 

where Pg±i/2 and Qg±i/2 are the radiative pressures and heat flux evaluated 
at the group interfaces. The Jacobian matrix J' of the hyperbolic system 
f fTSj) is given by 



P. 

J = -VuJ with J = ) ( =1 ^ - ^ ^ 7 1 (20) 




where ' denotes derivatives with respect to /. It can be shown that the 
trace and the determinant of J are both strictly positive. The eigenvalues 
of system (JT8l) are thus always of the same sign (i.e. opposite to that of 
P), which enables us to use a standard upwind scheme with respect to T) to 
calculate the values for P and Q at the group interfaces. This yields 



^9-1/2 



^g/AVg if P > 

X<,_i/Az/g_i if I? < 

(21) 

X,+i/Az/,+i ifl?>0 

^g/AVg if I? < 



where X = P or Q. This shows that the radiative energy and flux are advected 
from one group to the other depending on the sign of the velocity divergence. 
It is straightforward to show that the inclusion of these terms preserves the 
flux limitation condition |/| < 1 as long as < which is always true (see 
Fig.[T]). It would of course be possible to use higher order schemes to evaluate 
the quantities at the group interfaces by computing slopes using the usual 
methods. For the sake of conciseness this will not be explicited here. 



3. Numerical method and tests 

3.1. Numerical method 

In this section, we briefly present our global strategy to integrate the 
coupled RHD system (I7j)-(IHD (the method is identical to the one reported in 
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Gonzalez et al. [7|] apart from the terms involving the frequency differentials 
which were not included). In order to have a tractable time step, the radiative 
transport needs to be treated implicitly. However, it is most of the time more 
efficient to retain an explicit scheme for the hydrodynamics. We therefore 
use the following splitting scheme. 

In the first step the hydrodynamics system ([7]) is solved explicitly without 
the source terms. It is integrated using a classical second order MUSCL- 
Hancock scheme. In the second step, the radiation and the coupling terms 
are solved implictly. The terms in system (|8]) involving frequency derivatives 
are discretized as presented above and using the velocity divergence from the 
hydrodynamic solver. The velocity coming from the hydrodynamic solver 
is also used to discretize the other terms involving velocity derivatives. For 
the hyperbolic radiative term (two first left hand-side terms of system ([8])) 
we use a HLLC solver with an asyrnptotic preserving correction in order to 
recover properly the diffusion limit |3[. System ([8]) is solved implicitly with 
the source terms of system ([7]) using a Raphson-Newton procedure. 

r dtp + V-(pu)" = 

Step W dt{pu) + V-(pu®u + pI)" = (22) 
[ dte + V- (u(e + p))" = 



Step 2 < 



dtEg + V ■ F^+^ + V ■ (uEj^+i + (Pg : Vu)"+i 
- (Vu : P,)"+i = c(ap,e,(T) - 



dtF, + c^V ■ P:J+i + V ■ (u ® Fj"+i + (F„ ■ Vu)"+i 



- (Vu : Qg)"+i 



dte = -Y^ (^ciapgOgiT) - asgEg) - ((Tf<;/c)u ■ F^) 

9=1 

dtipu) = {{oFg/c)Y 



n+l (23) 



where 



P„ 



^9+1/2 



(9^(z/Py)(iz/ and 



^9-1/2 



d,{uQ,)du . (24) 



'^9-1/2 



The discretization of the divergence of a quantity U in cell i is done 
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following 



V.i/ = ^^ (25) 

where the state U* i is the value of the quantity lA (hydrodynamic or radia- 

*+2 

live) at the interface i + \ (between cell i and i + solution to the Riemann 
problem with left and right states and U^^i, respectively. For a first 



order scheme U. , ^ = U and U i = Wj+i. For a second order scheme, the 

values of U are lineraly extrapolated to the interfaces using local gradients. 

It is of course possible to solve systems (C])-® without splitting using a 
fully implicity scheme. We have tried this in one dimension and did not find 
any significant differences with the splitting scheme presented above. 

Calculating ipi and ip2 for every grid cell at every timestep is computa- 
tionally demanding due to the presence of logarithm and power functions, 
and we have thus tabulated the functions using 100 points which are read 
in once by the code at the beginning of a run. A specific value of (fi or ip2 
is then found using a Hermitian cubic spline interpolation which is very fast 
and accurate; the errors between the interpolated and the true values are less 
than 0.01% throughout. 

In this section we validate the method for multigroup RHD using a series 
of tests. For this purpose the numerical scheme skeched above have been 
implemented in a one- dimensional Lagrangian hydrodynamic code. 

The boundary conditions are implemented using two ghost cells at the 
edges of the grid. These ghost cells are filled using various physical con- 
straints such as null gradient, refiexive boundary or user imposed conditions. 

We use a step by step progression in our test sequence in order to verify 
each aspect of the method with increasingly complex problems. We first 
make sure that the multigroup transfer model (no hydrodynamics included) 
is equivalent to a grey model if the opacities are independent of frequency. 
We then test the multigroup aspect of the method with frequency dependent 
opacities. We compare the results of these tests to a well-established kinetic 
model which solves the equation of radiative transfer ([1]) directly j3| • Thirdly, 
we investigate the coupling of the radiative transfer to the gas motion using 
'frozen hydrodynamics'. Finally, we perform full RHD tests. 
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3.2. Marshak waves 

3.2.1. Classical grey Marshak wave 

Our first test is to check that the multigroup model reduces to a grey 
Ml model for a gas with a frequency-independent opacity. We run a Mar- 
shak wave simulation, where the gas inside the grid is at rest with a uni- 
form density p = 10"^ g cm~^, temperature T = 300 K in equihbrium 
with the radiation and opacity n = 1000 cm^ g~^ independent of frequency, 
noting that a = up. The specific heat capacity of the gas is set so that 
pCy = 10~^ erg g cm~^ K~^. The planar grid extends from to 20 cm, 
using 500 cells. Boundary conditions: the radiative energy inside the left 
and right ghost cells is that of a black body at 1000 K and 300 K, respec- 
tively. The radiative flux inside the left and right ghost cells is zero. We ran 
two simulations; the first using a single frequency group from to oo (grey 
model) and the second using five frequency groups evenly spaced between 
u = — 1.5 X 10^^ s~^ plus a sixth group to cover the range 1.5 x 10^^ s~^ 
to oo. The results are shown in Fig. |5] (solid lines), compared to the kinetic 
model (dashed lines), at a time t = 1.36 x 10^^ s. 

The radiative temperature inside a particular group is defined by 



In the top panel, the curves from the multigroup simulation (and the ki- 
netic model) are plotted. The curves representing the gas and radiative 
temperatures for the mono- and multigroup simulations were virtually in- 
distinguishable and we show the percentage difference between them in the 
bottom panel. Note that the differences remain below 0.5% throughout. This 
shows that the multigroup scheme consistently reduces to a grey model for 
frequency independent opacities. 

The kinetic model solves the equation of transfer directly using in this 
case 100 spatial zones, 64 directions and 64 frequency bins (all the results 
from the kinetic model have been tested for resolution convergence). For 
a moment model, the total radiative temperature (summed over all groups; 




(26) 



and the total radiative temperature is 




(27) 
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' 2 4 6 8 10 12 

Distance x (cm) 

Figure 2: Top panel: Gas and radiative temperatures in the grey Marshak wave test for 
n{v) = 1000 cm^ g~^ at time t — 1.36 x 10~^ s. The sohd curves are from the muthgroup 
Ml model and the dashed curves represent the kinetic model. The red curve marked T 
is the gas temperature and the green curve marked is the total radiative temperature 
(summed over all groups). The other coloured curves marked 1 to 6 represent the radiative 
temperatures inside each group. Bottom panel: percentage difference between the Mi 
grey and multigroup models for the gas temperature (red) and the radiative temperature 
(green). 

bright green) and the gas temperature (red) are in excellent agreement with 
their kinetic counterparts, illustrating the validity of the Mi model for ra- 
diative transfer and proving that the multigroup model consistently reverts 
to a grey model in the case of frequency-independent opacities. 

3.2.2. Multigroup Marshak wave with frequency dependent opacities 

As a second step, we consider a frequency variable opacity in order to 
assess its effect on the Marshak wave test. The setup is identical to the grey 
test above, but the opacities in the groups 1 to 6 are (in cm^ g~^) 1000, 750, 
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500, 250, 10 and 10 respectively. We also used 50 extra cells with steadily 
increasing widths at the right end of the grid in order to ensure that the 
radiation in the low-opacity groups does not have time to reach the right 
edge of the grid. The first 500 zones are the same as above, but the total 
grid size is ~ 9 m. The results are shown in Fig. [3] (solid lines). The gas 




Figure 3: Same as in Fig. [2] but in the case of a frequency dependent opacity. 

and radiation temperatures T and are different from the ones in the first 
test. The radiation in the groups with weak opacities (notably groups 5 and 
6) has crossed the entire grid and has heated the gas at the right edge (the 
gas temperature at that point is now 330 K). The radiation in the groups 3 
and 4 has also travelled further than in the previous test but not as far as 
groups 5 and 6. We note that the radiative temperature of group 1 at the 
right edge is slightly higher than in the previous test (just above 300 K as 
opposed to 275 K). Since its opacity is unchanged, this shows that the gas 
has been heated by the radiation in the other groups and has re-radiated 
some of its energy into group 1. The curves from the kinetic model (using 
400 cells, 100 directions and 512 frequencies) are also plotted (dashed lines). 
There is an extremely good agreement between the multigroup and kinetic 
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gas temperatures. The total radiative temperatures differ somewhat more 
than in the previous test and this difference is due to larger discrepancies in 
the low opacity groups 5 and 6. and are very close to their kinetic 
counterparts at the left edge of the grid but then drop rapidly and stabilize 
to a lower value. This is a boundary condition effect explained by the fact 
that when differences between the left and right fluxes are large (which is the 
case at the domain boundaries since the flux in the ghost cells is set to zero) 
the Ml model becomes less accurate. A solution to this issue would be to 
consider an additional third moment equation or to solve two half equations 
(one for the flux travelling towards the left and the other towards the right) 
for the radiative flux [6]. 

3.2.3. Multigroup Marshak wave with frequency and temperature dependent 



In our third test we use frequency variable opacities which also vary with 
temperature. The setup is identical to the multigroup test above, but the 
opacities are set to 



where Tq = 300 K and K^g in the groups 1 to 6 are the same as in the previous 
test, namely (in cm^ g~^) 1000, 750, 500, 250, 10 and 10, respectively. The 
grid setup is identical to the previous test. The results are shown in Fig. ID 
This time, as the gas temperature T increases, the opacity also increases 
and the radiation is absorbed much more rapidly, the effect being the most 
noticeable for groups 1 to 4 where opacities are high. We see once again an 
excellent agreement between the multigroup Mi model and the kinetic model 
which used 512 frequencies, especially for the gas temperature. 

3.3. Radiation traversing a region with strong velocity variations 

The aim of this next test is to study energy exchange between groups due 
to Doppler effects when strong velocity variations are present in the fluid. 
We perform this test in vacuum (p = k = 0). Radation is cast from the left 
side into the computation volume, with a black-body spectrum at = 1000 
K and a unit reduced flux. The size of the box is L = 10 cm for 50 cells. 



opacities 




3/2 



(28) 
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Figure 4: Same as in Fig. |3]but in the case of K{iy) = ko(i')(T/To)^/^. 
The velocity is set to obey the following law 

if X < xq 

if Xq < X < Xi 

if xi < X < X2 (29) 

if X2 < X < Xs 

if X > Xs 

where A = 5 x 10^ cm s~^, / = 6 cm, xq = 2 cm, xi = 3.5 cm, X2 = 6.5 
cm, X3 = 8 cm (see Fig. [5]). We used 20 equally spaced frequency groups 
in the range — ?■ 2 x 10^"^ Hz, plus a last group to hold frequencies in the 
range 2 x 10^^ — t- 00. The radiative temperature at the boundaries is kept 
constant at 1000 K, and the radiative reduced flux is maintained at / = 1. 
The system is left to evolve until stationarity is reached. 

The difference in radiative energies E^, between the fixed {u = 0) and 
the moving {u = A) regions is shown in Fig. El The circles are the group 
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A sin^ 
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Figure 5: Gas velocity as a function of x. 

OS 




Frequency v (Hz) 



Figure 6: Difference in radiative energies between a stationary {u = 0) and a moving 
{u = A) black body as a function of frequency. The solid line is the analytical solution, 
the circles are the numerical solution. 
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average numerical solution. The solid line is the analytical solution, which is 
obtained by applying a Doppler shift in frequency to the spectrum 



v' = ^v(\- (30) 

where 




It is clearly visible in Fig. |6] that due to the frequency shift of the black 
body spectrum, the first three frequency groups have gained energy while the 
remaining groups have lost energy. The discrepency between the analytical 
and numerical solution (both averaged within frequency groups) is of the 
order of one percent throughout. We have also performed the test with 10 
and 40 groups which did not change the errors significantly. 

3.J^. Velocity gradient with frequency dependent opacities 

In this test, we are interested in studying the ability of the code to handle 
frequency-variable opacities and Doppler shifts in a flow with strong velocity 
gradients. The size of the box is 1.0 cm for 100 cells and is initially filled with 
a gas at T = 3 K in equilibrium with radiation and a velocity u = Vx. As a 
first step, we set V = 0. In this test, the hydrodynamics are frozen and the 
gas density is set to p = l/{Cx) g cm~^ where C = 10^ (see below). The gas 
opacity varies with frequency: k(z/) = 100 cm^ g~^ for z/ < 2 x 10^^ s~^ and 
fi:(i/) = 1 cm^ g~^ for z/ > 2 x 10^'^ s"^, with a smooth transition between the 
two regimes of width Au = 4.5 x 10^ s~^ (see Fig. [7]). We use 20 frequency 
groups to sample the opacities. We then inject from the left hand side a 
radiation with a Gaussian intensity profile with a FWHM measuring 2/3 of 
the width of the opacities transition region which comprised the same energy 
as a 1000 K black body radiation. 

The radiative temperatures inside the separate groups are shown in Fig. [8] 
(top), where only the relevant groups are presented. The radiation in the first 
11 groups is rapidly absorbed by the gas which has a high opacity at these 
frequencies, while in the higher frequency groups, the radiation propagates 
rapidly in a quasi-transparent medium. The presence of radiative energy in 
groups 13 and above shows that the gas has been heated by the incoming 
radiation and has re-radiated some of its energy. Since the heated gas radiates 
as a black body, the radiation fills all the groups which are very narrow 
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Frequency v (Hz) 

Figure 7: Gas opacity as a function of frequency in cm^ (solid line). Intensity of the 
injected radiation (normalised; dashed). The FWHM of the Gaussian radiative intensity 
profile measures 2/3 of the width of the opacities transition region. The vertical black 
dotted lines represent the frequency groups identified by their numbers. 

compared to the width of a Planck curve at T = 1000 K. As the opacity is 
weak in the high groups, the radiation there can propagate freely towards 
the right edge of the box. The lower groups 1 to 3 are also filled by the black 
body radiation but their radiation cannot escape due to the strong opacities. 

In order to study the effects of velocity gradients on the radiation trans- 
port, we ran a simulation with the velocity gradient P = 10^ s~^. In this 
case, V = C and a permanent regime is achieved. The hydrodynamics are 
still frozen, which is justified by the fact that even the fastest gas would only 
have time to move a very small distance (5 x 10~^ cm) compared to the box 
size (1 cm) over the simulation time of 5 x 10"^^ s. The results are shown 
in Fig. [8] (bottom). We note that this time, only the radiation in the first 
8 groups is absorbed and the radiation in groups 9 and above propagates 
freely. This shows that the radiation in the intermediate groups 9, 10 and 
11 (covering the opacities transition region) was initially slightly absorbed. 
Then, as the Doppler frequency shift increases (due to the increasing veloc- 
ity), the radiation moves to groups of higher frequencies where the opacity 
is much lower and the radiation is thus able to escape freely. 
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Distance x (cm) 



Figure 8: Group radiative temperature at time t = 4.8 x 10 s in the gradient test for a 
null velocity (top) and for a gas velocity which increases linearly with distance (bottom). 



3.5. Astrophysical radiative shock 

The next test is to make sure that the muhigroup radiative transfer model 
is correctly coupled to the gas hydrodynamics. We ran a 'grey' radiative 
shock simulation using exactly the same parameters as in Gonzalez et al. [3] ■ 
The gas inside the computational domain is initially at rest with a uniform 
density of p = 7.78 x 10"^" g cm~^, temperature T = 10 K in equilibrium 
with the radiation and opacity k, = 0.39 cm^ g~^. The size of the box is 
1.0 X 10^^ cm. We give the gas at the left boundary a velocity of 20 km s~^, 
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which generates the propagation of a radiative shock travelhng towards the 
right. We use 500 equally spaced spatial zones and 6 frequency groups (5 
groups evenly spaced between u = — 7 x 10^"^ s~^ and the last group holds 
frequencies from 7 x 10^^ s~^ to oo). The results at three different epochs are 
shown in Fig. [9] (solid lines). We have also run the radiative shock test with 
SiNERGHYiD using only a single group and the results are shown in Fig. M 
(dashed lines). The temperature profiles are virtually indistinguishable, as 
illustrated by the difference AT between the grey and multigroup curves 
which is plotted below. The largest AT is ~ 40 K for a peak temperature 
of 4000 K, i.e. only one percent. Some small differences are visible in the 
radiative precursor. This shows that the multigroup scheme is consistent 
with the grey model. 

We also note that the curves are identical to the ones in Gonzalez et al. 
0], which shows that the implicit code correctly solves the equations of RHD. 




o 
o 




I 2xlo'° 4x10^° 6xlo'° 8xlo'° lo" 



Distance x (cm) 

Figure 9: Gas temperature in the radiative shock test as a function of distance at times 
t = 4.0 X lO'^ s (red), 7.5 x 10^ s (green) and 1.3 x lO'* s (blue) for the monogroup run 
(dashed) and the multigroup model using 6 groups (solid). The difference between the 
mono- and multigroup curves is shown in the bottom panel. 
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3.6. Multigroup radiative shock in xenon gas 

Our final test is to link sinerghyid to the ODALISC0 database of gas 
opacities in order to realistically model the evolution of a radiative shock in 
a xenon ^54Xe gas. The ODALISC database aims to provide spectral opacities 
as well as mean opacities (Rosseland and Planck) of many elements for a 
wide range of physical conditions. 

The gas inside the box is initially at rest with a uniform density of p = 
10"'^ g cm~^, temperature T = 1 eV in equilibrium with the radiation. The 
size of the box is 36 cm with 550 zones; the first 100 cells have logarithmically 
increasing sizes, the first 500 zones cover the range — 2 cm and the last 
50 have steadily increasing sizes, covering the range 2 — 36 cm. We give the 
gas at the left boundary a velocity of 60 km s~^, which generats a radiative 
shock travelling towards the right. We use an ideal gas equation of state with 
atomic mass number 131. 

The opacities for the Xe gas were taken from the ODALisc database 
(gomme average atom model). They depend on the gas temperature and 
density (often more strongly on temperature) as well as on the frequency. 
The opacity /t(z/) for the Xe gas for a density p = 10~^ g cm~^ and tem- 
perature T = 1 eV is shown in Fig. [10] (black solid curve). We used five 
groups to sample the opacities from u = 10~^ to 770 eV; the colour bands 
in Fig. [TU] illustrate the group decomposition of the frequency domain. Fre- 
quencies below 10^^ eV and above 770 eV are ignored, as gas temperatures 
in the box remain under 30 eV (except the very narrow temperature spike). 
The gas at such temperatures does not radiate strongly at these frequencies. 
We then computed the Planck (upg) and Rosseland (uRg) mean opacities in 
each group. 

Since the gas temperature and density evolve in time, the opacities need 
to be calculated at each timestep in each grid cell. The method we used 
to compute the opacities is to read in from the database a grid of opacities 
for the temperature range 0.01 to 100 eV and the density range 10~^ to 
0.3 g cm~^ at the start of the run. From this, we then compute Upg and Kp^g 
in each group at each point (p, T) which are stored into an array. During the 
simulation, a particular group opacity at any T and p is then found using a 
simple four-point interpolation using the array data. 

The gas and radiative temperatures for our simulation of a multigroup 



^ http: / /irfu. cea.fr/Projets / Odalisc / 
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Frequency v (eV) 



Figure 10: Xenon opacities at p = 10 g cm ^ and T — 1 eV as a function of frequencies. 
The colours illustrate the decomposition of the frequency domain into five groups from 
= IQ-^ to 770 eV. 

radiative shock in a Xe gas at a time t = 10~^ s are shown in Fig. [11] (top 
panel), along with the temperatures from an identical but grey run where 
only a single group over the same frequency range is used. A characteristic 
peak in the gas temperature (bright red) can be seen just around x = 6 x 10~^ 
cm at the shock. There is a strong radiative precursor (bright green) which 
extends all the way to x = 2 cm. The contributions to the precursor are 
clearly visible; at first the energy from group 5 contributes the most but 
subsequently gets dominated by group 4, 3 then 2 as we move further away 
from the shock. The radiation in the low-frequency group 1 does not appear 
to contribute to the dynamics of the shock. We note that the radiation in 
all the groups apart from group 2 gets absorbed fairly rapidly (none get 
past X = 0.6 cm), whereas since the opacities in group 2 are the lowest (see 
Fig. [TU]) . the radiation there propagates to a greater distance. Differences in 
the positions of the tip of the radiative precursors in the other groups further 
illustrate the effect of variable gas opacities. 

Let us note major differences between the mutigroup and the grey mod- 
els. Due to the fact that opacities are averaged over the entire frequency 
range (the low opacities are biased towards a higher value and vice versa), 
the radiation in the grey run suffers greater absorption far away from the 
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Figure 11: Top panel: Gas temperature T (bright red), total radiative temperature Tj. 
(summed over all groups, bright green) and individual group radiative temperatures (la- 
belled 1 to 5) in the multigroup simulation of a radiative shock in Xe gas at a time 
t = 10^^ s. We have also included the gas temperature (light blue) and radiative temper- 
ature (black) from the grey run. Bottom panel: Gas opacities for each group as a function 
of distance. The black curve represents the opacity in the grey run (averaged over all 
frequencies). 
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shock and its radiative precursor (black curve) does not extend as far as in 
the multigroup case. However, between 0.3 and 0.7 cm, the grey radiative 
temperature is much higher than the multigroup one. For this reason, the 
gas is heated to a greater extent and the grey gas temperature (light blue) is 
higher than the multigroup one around 0.3 — 0.4 cm. We also note that the 
multigroup gas temperature in the range 6 x 10^^ — 0.2 cm is higher than for 
the grey run. 

The opacities of the gas in each group and the opacity for the grey run 
are plotted as a function of distance in Fig. [11] (bottom panel). This is an 
excellent illustration of how the opacities are affected by the gas temperature. 
For instance, we see that the opacity in the first group is of the order of 
10^ cm^ g~^ in the cold gas ahead of the radiative precursor (right hand side) 
whereas it gains over three orders of magnitude in the hot post-shock gas. 
For groups 3, 4 and 5, the opposite occurs; the opacity is high before the 
shock and low after. The opacities in the temperature transition region has 
diverse behaviours. Most importantly, the curves are very different from the 
grey opacity (black curve) which is very constant with a single peak around 
0.4 cm corresponding to the jump in gas temperature (see light blue curve in 
the top panel). Interestingly, the pre- and post-shock grey opacities are very 
similar. It becomes very clear that the grey model cannot correctly represent 
the varied spectrum of opacities in such a situation, which are crucial to the 
evolution and dynamics of the shock. 

4. Conclusions 

We have developed a multigroup model for RHD using the Mi moment 
model. The equations of radiative transfer are solved in the comoving frame. 
In order to account for the opacity variations as a function of frequency, we 
introduced frequency groups and applied the Mi closure inside each of them. 
This gave rise to new terms depending on the frequency when coupled to the 
hydrodynamics. We use a finite volume method in the frequency domain in 
order to evaluate these new coupling terms which account for energy exchange 
between neighbouring groups due to the Doppler effect when strong velocity 
gradients are present in the gas flow. 

We have verified our method using a series of tests for both radiative 
transfer alone and radiative transfer coupled to hydrodynamics. In the case 
of the radiative transfer tests, the method was found to be successful in 
reproducing the results obtained with a kinetic code, at a much lower com- 



25 



putational cost. We have shown that the model reverts to a grey model for 
frequency independent opacities and that the model is capable of treating 
the effects of strong velocity gradients in a gas with frequency dependent 
opacities. 

Finally, we have coupled the sinerghyid code to the opacities from the 
ODALISC database to realistically simulate the propagation of a radiative 
shock in a Xe gas. We noted major differences between the multigroup and 
the grey models, showing the importance of accounting for the frequency 
variability of the gas opacities. The next step in this study will be to use 
a realistic equation of state for the Xe gas using the ODALisc database for 
more realistic simulations. ODALisc also has a number of different methods 
to calculate the opacities for each element. We will study the influence of the 
uncertainties of the opacities on the results of simulations of radiative shock 
in Xe in a following paper. In this work will also be included an investigation 
of the impact of the choice of frequency group boundaries and the number 
of groups on the results. 

We have now begun the implementation of this multigroup method for 
radiative transfer in the 3D radiation magnetohydrodynamics code HERA- 
CLES [7]. The development of such a tool for hydrodynamical simulations 
will prove extremely important for future studies, in particular, in the field 
of astrophysics where high-energy radiation from bright stars or supernovae 
gets absorbed by dense clouds which re-radiate the energy in the infrared, or 
inside dense stellar atmospheres where many chemical elements are present. 
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